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Abstract 

We present a grand canonical Monte Carlo simulation study of the phase diagram of a Lennard- Jones fluid adsorbed in 
a fractal and highly porous aerogel. The gel environment is generated from an off-lattice diffusion limited cluster-cluster 
aggregation process. Simulations have been performed with the multicanonical ensemble sampling technique. The biased 
sampling function has been obtained by histogram reweighting calculations. Comparing the confined and the bulk system 
liquid-vapor coexistence curves we observe a decrease of both the critical temperature and density in qualitative agreement 
with experiments and other Monte Carlo studies on Lennard- Jones fluids confined in random matrices of spheres. At variance 
with these numerical studies we do not observe upon confinement a peak on the liquid side of the coexistence curve associated 
with a liquid-liquid phase coexistence. In our case only a shouldering of the coexistence curve appears upon confinement. This 
shoulder can be associated with high density fluctuations in the liquid phase. The coexisting vapor and liquid phases in our 
system show a high degree of spatial disorder and inhomogeneity. 

PACS numbers: 61.20.Ja, 64.70.Fx, 05.70.Fh 



I. INTRODUCTION 

The phase behavior of fluids and fluid mixtures con- 
fined in porous and disordered materials represents a field 
of continuing theoretical and experimental interest due to 
a variety of applications in industrial technology 1 .. For 
example, porous materials are employed as adsorbents in 
many industrial processes, such as catalysis, adsorption 
separation, filtration and purification. 

Experimental studies on the phase behavior of fluids 
adsorbed both in highly porous materials, such as silica 
aerogel o 2 i 3 i 4 i 5 i 6 i 7 , and in less porous ones, such as Vy- 
cor glasses&2ii2iiiii2ii£ , have been performed. All these 
studies have shown that the phase diagram of the flu- 
ids confined in the porous structures is strongly altered 
in comparison with the corresponding bulk systems. For 
4 He& and NjL confined in very dilute aerogels it has been 
found that the liquid-vapor coexistence curve is much 
narrower and the critical temperature is lower than in 
the bulk. It has also been observed a shift of the criti- 
cal density towards the liquid phase and an increase of 
the coexistence vapor phase densities attributed to the 
attractive fluid-gel interactions. 

Very few theoretical studies on fluids confined in dis- 
ordered porous materials have been performed. A theo- 
retical approach based on considering the gel as a ran- 
dom field acting on the fluid succeeded in reproducing 
the important feature of spatial inhomogeneity- 14 . Im- 
portant phenomena like wetting cannot nonetheless be 
reproduced by such model. A model of "quenched- 
annealed" binary mixture was first studied by Madden 
and Glandt 15 and then widely used in dealing with the 
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problem of fluids phase separation in porous materials. 
Following this model, several integral equation theories 
and computer simulation studies have been successively 
proposed, in which the porous material is described by a 
random matrix of spheresiS*i2*i2ii2*2S. 

In particular the simulation work done by Page and 
Monson 19 dealt with the calculation of the phase diagram 
for a system representative of methane adsorbed in a sil- 
ica xerogel and that of Alvarez and aliSS has investigated 
the sensitivity of the confined fluid phase behavior to the 
matrix realization. The results of these studies about 
the vapor-liquid coexistence curve arc in good qualita- 
tive agreement with experiments. The novelty of these 
works not observed in experiments is the occurrence of 
an additional phase transition between a medium den- 
sity liquid and an high density one. The liquid-liquid 
phase coexistence has been associated with the wetting 
properties of the fluid in the more dense regions of the ad- 
sorbent. This second phase transition however has been 
found to be very sensitive to the matrix realization, while 
the liquid-vapor coexistence properties are seen to be ro- 
bust to variations in the solid structur o 20 ! 21 . This be- 
havior has been explained observing that the additional 
liquid-liquid transition is associated with the filling of low 
porosity regions in the matrix, which can occur for some 
configurations but not for others. The picture which 
emerges from all these calculations is that the coexis- 
tence vapor and liquid phases for the confined fluids are 
disordered and inhomogcneous as a consequence of the 
adsorbent structure randomness. 

The description of the aerogel as a random matrix of 
spheres however is not sufficiently realistic. The gel is or- 
ganized into a fractally correlated structure experimen- 
tally identified by a strong diffraction peak at small wave 
vectors measured by small angle x-ray or neutron scat- 
tering22i22i22i2£. The fractal behavior of aerogels is as- 
sociated with the irreversibility of the gel formation dy- 



namics, which proceed by random colloidal aggregation 
of silica particles. The fractal dimension of silica aerogels 
is about 1.80 in three dimensions 2 ^. 

Structures with similar fractal dimension have been 
generated by Monte Carlo simulations using several hi- 
erarchical cluster-cluster algorithms 2 ^' 2 ^. In particular 
the Diffusion Limited Cluster-Cluster Aggregation proce- 
dure (DLCA) is a gel growth process widely used in the- 
oretical studies of aerogels 2 **^. It has been shown that 
the structural properties and fractal dimension of DLCA 
gels well agree with experimental data on silica aerogels. 
A theoretical study of thermodynamical and structural 
properties of a Lennard-Jones fluid adsorbed in a highly 
porous DLCA aerogel by means of an integral equation 
approach has been recently proposed^. The phase dia- 
gram as well as the structural correlations of the confined 
fluid are found to be influenced by the specific properties 
of the gel, such as its connectivity and fractal behavior, 
particularly at low fluid densities. These results support 
the idea that a realistic modeling of the gel environment 
is necessary to deeply investigate the effects of the porous 
medium structural properties on the phase behavior of 
the adsorbate, especially on the second disorder-induced 
liquid-liquid transition. 

In this paper we present the results of a grand canoni- 
cal Monte Carlo simulation study about the phase behav- 
ior of a Lennard-Jones fluid confined in a dilute DLCA 
aerogel, where fluid and gel particles have hard core di- 
ameters of equal size and interact by means of an hard 
sphere potential. In order to locate the liquid-vapor co- 
existence curve we have calculated the density distribu- 
tion functions. The multicanonical ensemble sampling 
(MES) procedure has been employed to investigate ther- 
modynamical states in the subcritical region 28-29 30 . In 
the next section we describe the algorithm to build our 
confining system. In Sec. 3 we give details of our simula- 
tion. In Sec. 4 we present the results about the confined 
fluid phase diagram and compare it with that of the bulk 
system. Sec. 5 is devoted to conclusions. 



II. CONFINING SYSTEM 

The algorithm we used to generate the aerogel config- 
uration is the three dimensional off-lattice extension by 
Hasmy and al' of the DLCA procedure first proposed 
by Kolb and Hermann 2 !. The DLCA process is an itera- 
tive method which starts with a collection of N identical 
spherical particles of diameter a a randomly placed in a 
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cubic box of side L, with volume fraction rj = -jrc. 
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Aggregation proceeds via a diffusion motion of the par- 
ticles. 

If during their motion two clusters collide they stick 
together forming a new single aggregate. The process is 
terminated when a single cluster is obtained. Periodic 
boundary conditions are used at the edges of the simula- 
tion box. Details of the DLCA algorithm can be found 
in references 2 ^!. We have generated with this algorithm 
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FIG. 1: On the top plot of g(r) versus r/a a for an aerogel 
containing 515 particles in a box of edge L = 15<r a . This curve 
results from an average of 50 simulations. On the bottom 
Log-Log plot of S(q) versus qo a . 



configurations for aerogels containing 515 particles in a 
box of edge L = 15cr a , corresponding to a volume fraction 
i] = 0.08 and a porosity P, defined as the ratio between 
the free space volume and the volume occupied by the 
gel, of 92%. 

We have calculated for these aggregates the radial dis- 
tribution function g(r) and the static structure factor 
S(q) following the lines of the simulation work by Hasmy 
and aliS^ Our results are reported in Fig. ^ The ra- 
dial distribution function g(r) presents a strong peak at 
r = a a , associated with bonds between contacting parti- 
cles, and other features corresponding to short range cor- 
relations between particles belonging to the same cluster 
(see references' 2 ^ for further details). Before reaching 
the asymptotic limit of 1, g(r) exhibits a minimum cor- 
responding to distances between particles located at the 
boundary of the clusters. The location of this minimum 
gives an estimate of the mean clusters size £. In our case 
the mean cluster size is about 4cr a . As a consequence of 
the presence of the minimum in g(r), the scattering func- 
tion S(q), which is related to the Fourier transform of 
g(r) — 1, exhibits a pronounced peak at small q. The gel 
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exhibits a fractal structure in the intermediate q-range. 
The fractal dimension can be estimated from the power 
law behavior of the scattering function, S(q) ~ q~ D in 
this q-range. In our case we have obtained D ~ 1.74 
consistent with the value expected for silica aerogels in 
three dimensions 24 25 . 



III. SIMULATION DETAILS 

In our simulations the interactions between fluid par- 
ticles are described by a Lennard-Jones potential with a 
collision diameter a and a depth e truncated at r c = 2.5a. 
The aerogel particles are quenched and interact with the 
fluid particles by means of a purely repulsive hard sphere 
potential. We assume that fluid particles have a colli- 
sion diameter of the same size as the aerogel hard core 
a = a a . The edge of our simulation box is L = 15er. 
In the following Lennard-Jones units will be used: a for 
lengths, e for energies and e/fce for temperatures. 

In order to study the phase diagram of the Lennard- 
Jones fluid we have carried out Monte Carlo simulations 
in the grand canonical ensemble employing the algo- 
rithm used by Wilding 2 ^. MC moves consist of either 
an insertion or a deletion attempt, proposed with equal 
probability. Particles movements are implicitly imple- 
mented through insertion and deletion moves. 

In order to calculate the liquid and vapor coexistence 
densities the chemical potential can be varied, at fixed 
temperature, until a bimodal shape of the particles num- 
ber distribution P(N) is obtained. The coexistence den- 
sities at the selected temperature correspond to the two 
peaks positions. To obtain also the location of the liquid- 
gas saturation line in the (a*,T) plane the equal peak 
weight criterion 3 ^ can be applied. According to this 
criterion the chemical potential is tuned, at constant 
temperature, until the measured P(N) is double-peaked 
with equal area under the two peaks. The corresponding 
chemical potential belongs to the fluid coexistence curve. 

This conventional grand canonical ensemble Monte 
Carlo technique becomes however impractical in the sub- 
critical region we are interested in, due to the large free 
energy barrier separating the two phases and hindering 
spontaneous fluctuations of the system from the liquid 
to the gas phase and viceversa. In order to enhance 
crossing of this free energy barrier we have used the re- 
cently proposed MES 28,29-30 . With this powerful tech- 
nique sampling is made from a non Boltzmann distribu- 
tion, with a modified Hamiltonian H = H+g(N), where 
H(r, N) — E(y) — fiN is the configurational Hamilto- 
nian. The biased sampling function g(N) has to be cho- 
sen in such a way that the measured distribution P (TV) 
is nearly flat. In this way the mixed-phase configurations 
will be sampled with approximately the same probability 
as the gas and liquid configurations. The best choice for 
the biased sampling function would be g(N) — InP(N), 
where P(N) is the distribution we are looking for. An 
estimate P(N) can be obtained by extrapolation mak- 



ing use of the histogram reweighting technique 33 . In this 
approach we can accumulate the joint probability dis- 
tribution of system energy and particles number for a 
thermodynamical state characterized by a temperature 
To and a chemical potential /j,q (near the critical point 
where the interfacial tension is low)2S: 

e kT « V 

P(N,E\T , fM) ) = — (1) 

A> 

where V(N,E) is the density of states and Z (T ,/jIq) 
is the grand partition function. An estimate P of this 
distribution for another thermodynamical state (Ti,//i) 
can now be provided: 

Pi^ElTu^) = ^L e - { ^-^ ] P{N,E\T ^ ) (2) 

The new thermodynamical point (Ti, /ii) has to be suf- 
ficiently close to (To, /xo) m such a way that the statistical 
weight of the new configuration is not too different from 
the previous one. Integrating over the system energy E, 
we can finally obtain a suitable bias function : 

P(N\T 1 ,n 1 ) = J dEP(N,E\T u n{) (3) 

At the end of the Monte Carlo simulation at (Ti,/ii) 
with the biased Hamiltonian we obtain the almost flat 
P (TV). The real particles number distribution P(N) per- 
taining to (TifHi) is recovered from: P(N) = P (TV) ■ 
P(N). 

This method has given more accurate results than the 
well known Gibbs ensemble technique and Gibbs-Duhem 
integration scheme both for bulk and confined fluids 20 - 29 . 

We have performed extensive MES simulations in order 
to calculate the confined fluid phase diagram for a range 
of subcritical temperatures from T = 0.96 to T — 0.79. 
Our simulations involved about 200 ■ 10 6 steps for system 
equilibration and, depending on temperature, from 1-10 9 
to 2 • 10 9 steps for calculating the ensemble averages. 

Runs of this order of magnitude are necessary in order 
to supply enough indipendent samples from each of the 
two phases and consequently obtain density distribution 
functions with enough accuracy for our analysis 30 . We 
note that on decreasing temperature the rate of the in- 
terchanges between the liquid and gas phases becomes 
more rare and longer runs are needed. 

IV. RESULTS 

A. Simulation in the two-phase region 

First of all we need to locate the coexistence curve of 
the confined fluid and in particular the near-critical re- 
gion. For this sake we have calculated with standard 
grand canonical simulations adsorption isotherms of the 
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density p versus the chemical potential /j, for tempera- 
tures ranging from T = 0.80 to T = 0.95. Our results 
are shown in Fig. [3 FIGURA2 For a given temperature, 
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FIG. 2: Adsorption isotherms for a Lennard- Jones fluid 
in a DLCA aerogel for a range of subcritical temperatures 
(T = 0.80 diamonds, T = 0.85 triangles up, T = 0.90 circles, 
T = 0.92 triangles left and T = 0.95 squares). The fluid den- 
sity has been normalized by the void fraction (1 — 77) of the 
simulation box volume, where r\ is the gel volume fraction. 
All quantities are in Lennard- Jones units. 

we observe evidence of a coexistence at the same chemi- 
cal potential between a dilute gas and a medium density 
liquid in our confined system. We have not observed in 
our isotherms the occurrence of a liquid-liquid phase co- 
existence, at variance with previous computer simulation 
studies on the phase behavior of Lennard-Joncs fluids in 
disordered non fractal matrices of spheres 19,20 . 

In our calculations we have not considered desorption 
isotherms and the hysteresis associated usually with cap- 
illary condensation^^ since we are interested in the co- 
existence between equilibrium phases. 

We note that due to the slope of the liquid portion of 
the isotherms a small change of the chemical potential 
induces a strong variation in the liquid density. There- 
fore it is difficult to tune the chemical potential with the 
appropriate resolution to locate the coexisting densities 
with enough accuracy. On approaching the critical region 
this problem becomes even worst. 



B. Multicanonical ensemble sampling 

In order to better investigate the coexistence region 
and in particular the possible existence of a liquid-liquid 
phase coexistence we carried out extensive MES calcula- 
tions. 

From our isotherms calculations we have roughly esti- 
mated the region where the onset of the bimodal shape of 



the particles number distribution P(N) can be observed. 
In our analysis we started from the double-peaked P(N) 
obtained at the temperature T = 0.96. The histogram 
reweighting technique allowed us to estimate the coexis- 
tence liquid and vapor densities and chemical potential at 
this temperature. Then by making use of the MES pro- 
cedure we investigated a wide range of subcritical tem- 
peratures. In Fig.|31we report the distributions P(N) for 
the thermodynamical states investigated starting from 
T = 0.96. 
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FIG. 3: Coexistence particles number distributions for sub- 
critical temperatures ranging from T — 0.96 to T = 0.79 in 
Lennard- Jones units. Curves with closer peaks correspond to 
higher temperatures. 

At higher temperatures the liquid and gas have close 
densities and the double-peaked P(N) distribution shows 
two partially overlapped peaks. As the temperature de- 
creases the difference between the coexisting vapor and 
liquid densities becomes more marked. 

In the intermediate temperature regime we do not ob- 
serve the occurrence of a third peak associated with a 
second coexisting liquid phase. We will comment more 
extensively about this point in the following. Table |U 
summarizes the temperatures at which we performed our 
simulations together with the values of the coexistence 
densities. The MES algorithm allowed us to follow the 
system down to T — 0.79, where the liquid density is 
more than f times higher than the vapor one. Snapshots 
from our simulations for the coexistence liquid and va- 
por configurations at one temperature are shown in Fig.|3] 
and Fig. In both the coexisting phases the spatial dis- 
tribution of the fluid molecules is inhomogeneous and the 
pores of the aerogel are far from being uniformly filled by 
the fluid. 

Fig. shows the temperature versus density phase di- 
agram we obtained from the peaks locations in the P(N) 
distributions depicted in Fig. |3| compared with that of 
the bulk fluid 29 . We find that the confined fluid phase 
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FIG. 4: A snapshot for the gas phase near the coexistence at 
T = 0.915. The light grey and the black spheres represent the 
gel and fluid particles, respectively. The boxlength is L — 15. 
All quantities are in Lennard- Jones units. 



diagram is substantially modified by the presence of the 
aerogel : both the critical temperature and density are 
lower than in the bulk and the range of the vapor-liquid 
coexistence curve is much less extended. We recall that 
in the bulk the critical parameters are: T c — 1.1876, 
Pc = 0.3197 and [i c = -2.7782a. 

Our findings about the gas-liquid coexistence prop- 
erties are in qualitative agreement with previous com- 
puter simulation studies on Lennard- Jones fluids con- 
fined in random spheres matrices with purely repulsive 
adsorbent-adsorbate interactionsii2i2ii. 



FIG. 5: A snapshot for the liquid phase near the coexistence 
at T = 0.915. The light grey and the black spheres represent 
the gel and fluid particles, respectively. The boxlength is 
L — 15. All quantities are in Lennard- Jones units. 
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FIG. 6: Phase diagram for the confined fluid obtained from 
the particle number distribution peak locations (circles) com- 
pared with that of the bulk fluid-S (stars). All quantities are 
in Lennard- Jones units. 

None-the-less the phase diagram reported in Fig. 
does not present a liquid-liquid coexistence region. It 
shows instead a well defined shoulder on the liquid side 
boundary for intermediate temperatures in the range in- 
vestigated. 

We must however stress that confining primary par- 
ticles in real aerogels are substantially larger than the 
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FIG. 7: Blow up of the liquid peaks of the number particles 
distributions depicted in Fig. The bold solid lines corre- 
spond to the liquid peaks generating the observed shoulder 
in the phase diagram (T = 0.885, T = 0.87, T = 0.86 and 
T = 0.85 in Lennard- Jones units). 



fluid molecules, therefore a study on a much larger sys- 
tem could help to shed light on this issue. 

In order to better define the nature of the liquid 
phase for the thermodynamic points corresponding to 
this shouldering, in FigQwe show a blow up of the liquid 
peaks of the P(N) distributions centered around the in- 
termediate temperature region. The liquid peaks gener- 
ating the shoulder in the coexistence curve are indicated 
with bold solid lines. The liquid peaks shape shows a 
peculiar behavior with varying temperature. For tem- 
peratures ranging from near the critical point to about 
T = 0.90, the shape of the liquid peaks is symmetric 
around the peaks maximum as observed for the bulk 29 . 
For the reduced temperature corresponding to the on- 
set of the shoulder, T — 0.885 (the first bold curve on 
the left), the liquid peak shows a bump on its right side 
that renders its shape slightly asymmetric. We observe 
that the position of this bump roughly corresponds to the 
liquid peak maximum of the next P(N) distribution, cor- 
responding to T = 0.87. For T = 0.87 (the second bold 
curve starting from the left) it is conversely observed a 
bump located on the left side of the liquid peak at nearly 
the same position of the peak maximum measured at 
T = 0.885. We additionally note that the main peaks 
of these two curves show a more pronounced separation 
from each other compared with the other curves, in spite 
of being the temperature jump between one curve and 
the next approximatively similar for all the curves inves- 
tigated. The next two liquid peaks , corresponding to 
the end of the shouldering, T = 0.86 and T = 0.85, are 
less asymmetric and closer than the two previous ones. 
For temperatures lower than T — 0.85, after the shoul- 
der in the phase diagram, the shape of the liquids peaks 



FIG. 8: The confined fluid liquid-gas saturation line in the 
(/Lt,T) plane for temperatures ranging from T = 0.96 to T — 
0.79. All quantities are in Lennard- Jones units. 



returns to be symmetric and the liquid densities increase 
regularly on lowering temperature, similar to the high 
temperatures range. 

From these observations we can infer that at high tem- 
peratures the gas phase coexists with a medium density 
liquid. Correspondingly the coexistence curve decreases 
regularly with temperature. In the intermediate temper- 
ature range the liquid phase starts showing high density 
fluctuations and appears to coexists with a slightly more 
dense liquid. This region is characterized by an asymme- 
try in the liquid peaks of the P(N) and correspondingly a 
shouldering of the phase diagram, that mark a crossover 
to a slightly more dense phase. For lower temperatures 
the higher density liquid becomes the thermodynamically 
favoured phase coexisting with the vapor and the coex- 
istence curve resumes a regular descent. 

Finally we report the liquid-vapor coexistence curve in 
the (n,T) plane as obtained from our multicanonical sim- 
ulations. According to the equal peak weight criterion, 
described in Sec. 3, the location of the liquid-vapor sat- 
uration line in this plane can be obtained looking for the 
values of temperature and chemical potential at which 
the bimodal distributions P{N) have equal area under 
the two peaks. 

Fig. |H1 shows the temperature behavior of the coexis- 
tence chemical potentials at which we obtained the P(N) 
distributions reported in Fig. [31 Analogous to the bulk 
we observe a monotonic regular behavior. Similar to the 
temperature versus density phase diagram, there is a sub- 
stantial shrinkage of the range of chemical potentials at 
which we observe the coexistence of the liquid and gas 
phases with respect to the bulk. In fact in the bulk 2 ^ for a 
similar interval of temperatures as ours, i.e. AT = 0.17, 
a range A/i ~ 0.6 is found, while in our confined system 
A/i ~ 0.09, more than six times lower. 
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V. CONCLUSIONS 

We have performed a computer simulation study of 
the phase behavior of a Lennard- Jones fluid adsorbed in 
a highly porous fractal aerogel. The gel environment has 
been generated with the DLCA algorithm2£i2£ in order 
to obtain a more realistic confining structure than those 
reported in literature. Besides the confined fluid phase 
diagram has been calculated by performing multicanoni- 
cal ensemble simulations in the framework of the Monte 
Carlo grand canonical ensemble techniqueS&SSiSl. 

The proper bias has been found by applying the single- 
histogram reweighting technique^. This procedure is 
particularly suitable for analysing the transition state of 
the first order phase transitions in great detail. 

We found that the phase diagram of the confined fluid 
is substantially modified with respect to the bulk by the 
presence of the gel: the critical temperature and density 
are lower and the liquid-gas coexistence curve is much 
narrower than in the bulk. Our findings about the vapor- 
liquid coexistence properties are in qualitative agreement 
with experimental observations and previous theoretical 
studies. We have not found a clear evidence for a second 
liquid-liquid coexistence region in the fluid phase dia- 
gram. The presence of a liquid-liquid phase coexistence 
has been sometimes reported, depending on the partic- 
ular configuration generated for the adsorbent, random 
matrices of spheresi2i2£. Our coexistence curve shows a 
shoulder on the liquid side boundary in the intermedi- 
ate temperature range. For the same temperatures the 
liquid peaks shape of our particles number distributions 



are highly asymmetric and show a bump. However, since 
we do not find two distinctly resolved peaks, we must 
conclude that in our case coexistence between two liquid 
phases is never observed. We do observe instead a cross- 
ing of the liquid from a lower to an higher density phase 
upon decreasing temperature. 

We have employed in this study a very sophisticated 
technique for an accurate location of the coexistence 
curve. Due to the fractal structure of the confining envi- 
ronment the results should in principle not depend on the 
realization. However we cannot exclude that finite size 
effects could be present. Only extensive analysis with in- 
creasing sizes can confirm how sensitive the results are 
to the matrix realization. 

In order to gain deep insight about the modification 
of the phase diagram of fluids confined in aerogels the 
modeling of the adsorbent environment can be further 
improved to make contact with the experimental situa- 
tions. It would be interesting to explore in a more care- 
fully designed confining structure with a network of mul- 
tishaped voids if the two liquid phases will become dis- 
tinguishable. In particular the liquid-liquid coexistence is 
smeared out in our present simulation, where we observe 
only a shouldering in the phase diagram. It would be 
valuable to perform simulations on even larger systems 
to better investigate the role of the fractal behavior of 
the confining medium and perform an extensive analysis 
of finite size effects. 

Work on these improvements to our present study is 
currently in progress. 
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